Using Satellite Data for Species Distribution Modeling with GRASS GIS and R

Importing species records

GRASS code

# Import records
v.import input=aedes_albopictus.gpkg 
  output=aedes_albopictus

# List raster maps
g.list type=raster
r.colors map=lst_2014.150_avg color=celsius

# Display records
d.mon wx0
d.rast lst_2014.150_avg
d.vect aedes_albopictus icon=basic/circle \
  size=7 fill_color=black

Map

You can also get the occurrences directly from GBIF into GRASS by means of v.in.pygbif.

Creating random background points

GRASS code

# Create buffer around Aedes albopictus records
v.buffer input=aedes_albopictus output=aedes_buffer distance=2000

# Set computational region
g.region -p raster=lst_2014.001_avg

# Create a vector mask to limit background points
r.mapcalc expression="MASK = if(lst_2014.001_avg, 1, null())"
r.to.vect input=MASK output=vect_mask type=area

# Subtract buffers from vector mask
v.overlay ainput=vect_mask binput=aedes_buffer operator=xor output=mask_bg

# Generate random background points
v.random output=background_points npoints=1000 restrict=mask_bg seed=3749

Map

See extra slides for details about computational region and masks in GRASS GIS.

Create daily LST STRDS

GRASS code

# Create time series 
t.create type=strds temporaltype=absolute \
  output=lst_daily title="Average Daily LST" \
  description="Average daily LST in degree C - 2014-2018"

# Get list of maps 
g.list type=raster pattern="lst_201*" output=list_lst.csv

# Register maps in strds  
t.register -i input=lst_daily file=list_lst.csv \
  increment="1 days" start="2014-01-01"

# Get info about the strds
t.info input=lst_daily

Output

See t.create and t.register

Generate environmental variables from LST STRDS

Long term monthly avg, min and max LST

GRASS code

for i in $(seq -w 1 12) ; do 
  t.rast.series input=lst_daily method=average \
    where="strftime('%m', start_time)='${i}'" \
    output=lst_average_${i}
  t.rast.series input=lst_daily method=minimum \
    where="strftime('%m', start_time)='${i}'" \
    output=lst_minimum_${i}
  t.rast.series input=lst_daily method=maximum \
    where="strftime('%m', start_time)='${i}'" \
    output=lst_maximum_${i}  
done

Output

See t.rast.series manual for further details.

Bioclimatic variables

GRASS code

# Install extension
g.extension extension=r.bioclim
 
# Estimate temperature related bioclimatic variables
r.bioclim \
  tmin=$(g.list type=raster pattern="lst_minimum_??" separator=",") \
  tmax=$(g.list type=raster pattern="lst_maximum_??" separator=",") \
  tavg=$(g.list type=raster pattern="lst_average_??" separator=",") \
  output=worldclim_ 

# List output maps
g.list type=raster pattern="worldclim*"

Output

Map

See r.bioclim manual for further details.

Spring warming

GRASS code

# Annual spring warming: slope(daily Tmean february-march-april)
t.rast.aggregate input=lst_daily output=annual_spring_warming \
  basename=spring_warming suffix=gran \
  method=slope granularity="1 years" \
  where="strftime('%m',start_time)='02' or \
         strftime('%m',start_time)='03' or \
         strftime('%m', start_time)='04'"

# Average spring warming
t.rast.series input=annual_spring_warming \
  output=avg_spring_warming \
  method=average

Map

See t.rast.aggregate manual.

Autumnal cooling

GRASS code

# Annual autumnal cooling: slope(daily Tmean august-september-october)
t.rast.aggregate input=lst_daily output=annual_autumnal_cooling \
  basename=autumnal_cooling suffix=gran \
  method=slope granularity="1 years" \
  where="strftime('%m',start_time)='08' or \
         strftime('%m',start_time)='09' or \
         strftime('%m', start_time)='10'"

# Average autumnal cooling
t.rast.series input=annual_autumnal_cooling \
  output=avg_autumnal_cooling \
  method=average

Map

Number of days with LSTmean >= 20 and <= 30

GRASS code

# Keep only pixels meeting the condition
t.rast.algebra -n \
  expression="tmean_higher20_lower30 = if(lst_daily >= 20.0 && lst_daily <= 30.0, 1, null())" \
  basename=tmean_higher20_lower30 suffix=gran nproc=7

# Count how many times per year the condition is met
t.rast.aggregate input=tmean_higher20_lower30 \
  output=count_tmean_higher20_lower30 \
  basename=tmean_higher20_lower30 suffix=gran \
  method=count granularity="1 years"

# Average number of days with LSTmean >= 20 and <= 30
t.rast.series input=count_tmean_higher20_lower30 \
  output=avg_count_tmean_higher20_lower30 method=average

Map

See t.rast.algebra manual for further details.

Number of consecutive days with LSTmean <= -2.0

# Create annual mask
t.rast.aggregate input=lst_daily output=annual_mask \
  basename=annual_mask suffix=gran \
  granularity="1 year" method=count

# Replace values by zero
t.rast.mapcalc input=annual_mask output=annual_mask_0 \
  expression="if(annual_mask, 0)" \
  basename=annual_mask_0

# Calculate consecutive days with LST <= -2.0
t.rast.algebra \
  expression="lower_m2_consec_days = annual_mask_0 {+,contains,l} \
  if(lst_daily <= -2.0 && lst_daily[-1] <= -2.0 || \
  lst_daily[1] <= -2.0 && lst_daily <= -2.0, 1, 0)" \
  basename=lower_m2_ suffix=gran nproc=7

# Inspect values
t.rast.list input=lower_m2_consec_days \
  columns=name,start_time,end_time,min,max

# Median number of consecutive days with LST <= -2
t.rast.series input=lower_m2_consec_days \
  output=median_lower_m2_consec_days method=median

We have all these maps in GRASS, how do we connect with R now?